home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / testdemo.demo < prev    next >
Text File  |  1997-07-08  |  7KB  |  203 lines
  1. ; $Id: testdemo.demo,v 1.4 1996/12/18 00:12:36 ali Exp $
  2.  
  3. ; IDL Demonstration script.  A sequence of IDL statements just
  4. ;    as they would be entered from the keyboard.
  5. ;
  6. if !d.n_colors gt 2 then loadct,5    ;Color display?
  7.  
  8. ;
  9. ;    Simple Vector demonstration:
  10. ;
  11. y = [0,1,2,3,4,5]    ;Define a 6 element integer vector
  12. plot,y,title='Simple Plot'  ;Plot it
  13. oplot,sqrt(y),color=151    ;overplot the square root
  14. ;
  15. n = 100            ;Size of vectors
  16. ;
  17. ;    n element vector demonstration:
  18. ;
  19. x = findgen(n)        ;n element floating array: 0, 1, 2, ..., n-1.
  20. plot,x,title='Original Vector' ;Show it
  21. y = sin(8*!pi/n * x) / exp(x/(n/2.))    ;Make a dampend sine wave
  22. plot,y,title='Dampened Sine Wave'
  23. for i=1,3 do oplot,y/(i+1),linestyle=i    ;Show differing line styles
  24. ;
  25. ;    Polygon filling:
  26. ;
  27. xx = [x, reverse(x)]    ;Vector of X values concatenated with reversed vector.
  28. yy = [y, -reverse(y)]    ;Same with Y but negate last half
  29. for i=0,3 do polyfill,xx,yy/(i+1),color=255/(i+1) & end
  30. ;
  31. ;    Example of spline interpolation.
  32. ;
  33. window,2,title='Splines'    ;make a new window
  34. a = randomn(aa,12)    ;12 element normal dist random vectors
  35. plot,a,title='Spline Interpolation' ;simple plot
  36. xx = findgen(121)/10.    ;grid for spline interpolation, 121 pts from 0 to 12.0
  37. yy = spline(findgen(12),a,xx) ;apply spline function to random points
  38. oplot,xx,yy,thick = 2    ;overplot with double thick lines
  39. ;
  40. ;    Gaussian fitting demo:
  41. ;
  42. window,1,title='Gaussian Fitting'
  43. x = findgen(50)        ;5 element vector
  44. v = randomu(aa,3)*[30,10,10]+[10,1,1] ;random #'s for center width and height.
  45. z = v[2] * (exp(-((x-v[0])/v[1])^2) + 0.2*randomn(aa,50)) ;make signal
  46. plot,z,title='Gaussian',psym=3
  47. zz = gaussfit(x,z,a)    ;Fit it using library function
  48. oplot,zz,color=200
  49. xyouts,2,.8*!y.crange[1],'!8y = a!I0!N e!A-z!E2!A/2!3',siz=2.5
  50. print,'Center',v[0],', 1/e width',v[1],', Height',v[2]
  51. print,'Calculated:'
  52. print,'Center',a[1],', 1/e width',a[2],', Height',a[0]
  53. ;
  54. ;
  55. ;    Frequency domain filter demonstration:
  56. ;
  57. window,2,xs=700, ys=500, title='Low Pass Filtering'
  58. !p.multi = [0,2,2]        ;Gang plots, 2 by 2.
  59. ynoise = y + randomn(aa,n)/4.    ;Add noise to corrupt dampened sine wave
  60. plot,ynoise, psym=3,title='Signal + Noise', color = 201
  61. oplot,y                ;Original signal as a line
  62. ;
  63. z = abs(fft(y,1))        ;Original power spectrum
  64. znoise = abs(fft(ynoise,1))    ;Corrupted Power spectrum
  65. ;                Log Plot of power spectra:
  66. ;                Plot original power spectum
  67. plot_io,shift(z,n/2),title='Power Spectra',ytitle='Log Power'
  68. ;                ;If color, fill original
  69. if demo$color then polyfill,findgen(n),shift(z,n/2)
  70. oplot, shift(znoise,n/2), psym=3, col = 201 ;Overplot noise
  71. ;
  72. freq = findgen(n)/n        ;Make a vector proportional to freq
  73. freq = freq < (1.-freq)        ;symmetrical about Nyquist freq
  74. filter = 1./(1 + (freq/(4./n))^2) ;Low pass butterworth
  75. plot, shift(filter,n/2), title='Butterworth Filter Function', color = 151
  76.  
  77. result = fft(filter * fft(ynoise, -1), 1) ;apply filter
  78. plot, ynoise, psym=3, title='Filtered Result',  color=201
  79. oplot, result, col=151        ;filtered result
  80. oplot, y, lines=2        ;original signal
  81.  
  82. !p.multi = 0            ;Reset standard plots
  83.  
  84. ;
  85. ;    Two dimensional array Section
  86. ;
  87. loadct,3            ;red color tables
  88. wset,0 & wshow,0        ;clean up things
  89. if !d.n_colors le 2 then top = 255 else top = !d.n_colors-1 ;Max Col index
  90.  
  91. a=shift(dist(40),20,20)        ;Make distance array
  92. a=exp(-(a/8.)^2)        ;Yes, the proverbial Gaussian.
  93.                 ;Demonstrate contour plot
  94. contour, /follow, a, title='Contour Plot',c_line=[0,1,2,3]
  95. threed,a,title='Stacked Row Plot'    ;Show stacked row plot
  96.  
  97. window, 1, title='Second IDL Window',xsiz=480,ysiz=360    ;make a new window
  98. surface, a, /save, /xst, /yst, bot=129, skirt=0. ;combine them
  99. contour, a, /t3d, /noerase, zval=1.0, /xst,  /yst, c_color=[101,151,201,251]
  100. ; More Demonstrations of Contouring:
  101. a = randomu(seed, 6, 5)        ;Make a 6x5 array of random numbers
  102.  
  103.  
  104. ;  Show some contours
  105. levels = findgen(9)/10. + .1
  106. contour, a, level=levels, title='Simple Contour Plot'
  107. b = min_curve_surf(a)
  108. wset, 0
  109. tek_color
  110. contour, b, level=levels, title='Minimum Curvature Surface'
  111. contour, b, /fill, level = levels, title='Polygon Fill', c_color=indgen(9)+2
  112. contour, b, color = 0, /overplot, /downhill, levels = levels
  113. wait, 1 
  114.  
  115. if not demo$color then demo$done = 1    ;bail out now if can't display images.
  116. ;
  117. ; Image processing demo.
  118. ;
  119. openr,1, demo$imagedir + 'cereb.dat',err=i  ;Open image file
  120. if i ne 0 then begin print,'Image files not found' &  demo$done=1 & end
  121. file = assoc(1, bytarr(512,512))    ;contains 512 by 512 byte images
  122. a = file[0]                ;read 1st image
  123. loadct, 3
  124. wset, 0                    ;use original image
  125. erase                    ;clean out display
  126. tv, a                    ;display image
  127. h = histogram(a)            ;get pixel distribution
  128. wset, 1                    ;other window
  129. plot, h[1:*], title='Pixel Distribution Histogram'
  130. wset, 0
  131. h_eq_ct, a            ;show histogram equalized image
  132. xyouts, 310, 480, /noclip, 'Histogram Equalized Image', /dev
  133. ;
  134. ;    Image subtraction.  The first image is before Iodine dye injection,
  135. ;        and the second image is after.
  136. ;
  137. b = bytscl(fix(a) - file[1],top=top)    ;subtract 2nd image
  138. close, 1                ;done with file
  139. stretch                    ;restore normal color table
  140. tv, b                    ;show difference image
  141. xyouts, 310, 480, /noclip, 'Mask Subtraction Image',  /dev
  142. b = hist_equal(b, top = top)        ;histogram equalize pixels
  143. tv,b                    ;Show it
  144. shades = rebin(b,64,64)            ;shades for later use (4d display)
  145. height = rebin(a,64,64)-50>0        ;Heights
  146. ;
  147. ;    Unsharp masking
  148. ;
  149. b = bytscl(fix(a) - smooth(a, 5),top=top) ;subtract smoothed orig from smoothed
  150. stretch                    ;restore color tables
  151. tv, b                    ;show diff
  152. xyouts, 330, 490, /noclip, 'Unsharp Masking', /dev
  153. h_eq_ct,b                ;Histogram equalize it
  154. ;
  155. ;    2D Fourier Filtering
  156. ;
  157. openr,1, demo$imagedir + 'abnorm.dat'    ;Open image file
  158. a=bytarr(64,64)            ;define array
  159. readu,1,a            ;read image
  160. close,1                ;done with file
  161. erase  & stretch        ;clean out display
  162. tvscl,rebin(a,256,256),0    ;show it
  163.  
  164. xyouts,0,480,/dev,/nocl,'Original Image'
  165. b = fft(a,1)            ;forward transform
  166. z = alog(abs(b))        ;Log power spectrum
  167. tv,rebin(shift(bytscl(z,top=top),32,32),256,256),1    ;show it
  168.  
  169. xyouts,256,480,/dev,/nocl,'Power Spectrum'
  170. freq = dist(64)            ;make a butterworth low pass filter
  171. filter = 1./(1 + (freq/5)^2)    ;cutoff = 5 cycles/image.
  172. tvscl,rebin(bytscl(shift(filter,32,32),top=top),256,256),2    ;show it
  173.  
  174. xyouts,0,224,/dev,/nocl,'Butterworth Filter Function'
  175. result = fft(filter * b,-1)    ;filter & retransform
  176. tv,rebin(bytscl(result,top=top),256,256),3 ;show result
  177. xyouts,256,224,/dev,/nocl,'Low Pass Filtered'
  178. ;
  179. ;    Combined view of image:
  180. ;
  181. window,2,title='Combined Display',xsize=400, ysize=400
  182. show3,smooth(a,5),/interp,sscale=2    ;show image 3 ways
  183. ;
  184. ;    Shaded view of image:
  185. ;
  186. window,1,title='Shaded Surfaces',xsize=480, ysize=400
  187. shade_surf,rebin(a,32,32)    ;Make a shaded surface
  188. ;
  189. ;    4D view of brain and circulation:
  190. ;
  191. ;    The elevation is the X-ray density, and the shading shows the 
  192. ;    blood perfusion.
  193. ;
  194. wset,0
  195. shade_surf,height, shades=shades, ax=70,/xst,/yst ;show it
  196. surface,rebin(height, 32,32),ax=70,col=0,/xst,/yst,/noer ;overprint lines
  197. loadct,15        ;Load a pretty color table
  198. xyouts,10,480,/dev,/nocl,font=0,'Elevation = X-Ray Density'
  199. xyouts,10,450,/dev,/nocl,font=0,'Shading = Blood perfusion'
  200. ;
  201. ;
  202. ;    ALL DONE!!
  203.